Microscopic model of Purcell enhancement in hyperbolic metamaterials 
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We study theoretically a dramatic enhancement of spontaneous emission in metamaterials with 
the hyperbolic dispersion modeled as a cubic lattice of anisotropic resonant dipoles. We analyze the 
dependence of the Purcell factor on the source position in the lattice unit cell and demonstrate that 
the optimal emitter position to achieve large Purcell factors and Lamb shifts are in the local field 
maxima. We show that the calculated Green function has a characteristic cross-like shape, spatially 
modulated due to structure discreteness. Our basic microscopic theory provides fundamental insights 
into the rapidly developing field of hyperbolic metamaterials. 
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I. INTRODUCTION 

Engineering light-matter interaction in nanostructured 
environment has recently been the focus of active stud- 
ies [IH3- in particular, the so-called hyperbolic meta- 
materials described by an uniaxial medium where the 
main components of dielectric and magnetic tensor have 
different signs, have attracted significant attention [5J- 
HOj . Realizations of the regime of the hyperbolic medium 
with negative components of the dielectric tensor have 
been reported for magnetized plasma [TT], graphite [12] , 
for metamaterials based on nanorod assemblies [13-15 
and for layered metal-dielectric structures [TBI H7j . In 
this regime, light wavevectors at a given frequency fill a 
surface of a hyperbolic shape, so that the area of hyper- 
bolic isofrequency surface, giving the photonic density of 
states, is infinite. As the result, the spontaneous emis- 
sion rate becomes infinite in the case of ideal hyperbolic 
medium 

Experimental reports on the Purcell factor enhance- 
ment, describing the spontaneous emission rate modi- 
fication, in hyperbolic metamaterials are already avail- 
able |U QI5HH]- Theoretical studies for various models 
have been also performed [52]. It has been shown that 
the Purcell factor should not actually diverge. It is rather 
determined by a cutoff in the wavevector space, stem- 
ming from spatial inhomogeneity of the medium |23H25| , 
a finite distance from the source to the medium [26] [27] , 
nonlocality of dielectric response [28] , or finite size of the 
emitter gj3]. 

The basic solid state model of either natural or artifi- 
cial material is a periodic lattice of unit cells. We adopt 
this model and consider hyperbolic material as the infi- 
nite cubic crystal of interacting resonant point dipoles. 
Similar models have been developed for various systems 
including lattices of quantum dots [30], optical atomic 
lattices [3TJ [32], 7-ray resonant nuclear scattering [33J 



as well as the discrete-dipole-approximation of the light 
scattering theory [34]. This general approach, despite 
certain approximations, has been successfully applied to 
the lattices of split-ring resonators [35] [35] . 

In this paper, we study optical properties of the infinite 
cubic crystal of resonant interacting point dipoles polariz- 
able only in one direction (see Fig. [IJ . This model allows 
us to reproduce the hyperbolic isofrequency surfaces of 
the uniaxial anisotropic metamaterials and accounts for 
the discrete character of metamaterials. Within this mi- 
croscopic model of hyperbolic metamaterial, we investi- 
gate the influence of emitter position within the unit cell 
of the metamaterial on its radiation properties. 

The paper is organized as follows. Section |Tl] outlines 
our theoretical model and approach. Calculated disper- 
sion and lattice Green function are discussed in Sec. Mil 
Section |IV| is devoted to the numerical and analytical 
results for the Purcell factor and Lamb shift in metama- 
terials with hyperbolic dispersion. 



II. DISCRETE DIPOLE MODEL 

We consider an infinite periodic cubic lattice rj of the 
point dipoles, characterized by the period a, and embed- 
ded in vacuum. Our approach can be straightforwardly 
generalized to allow for background dielectric constant 
e ^ 1. The emitter inside the lattice is modeled by a 
radiating dipole po is placed at the point tq. Structure 
geometry is sketched on Fig. [T] The self-consistent elec- 
tric field satisfies the following equation 

V x V x E - q 2 E = Anq 2 P , (1) 

where q — lo/c is the wavevector at the frequency u). The 
quantity P in Eq. (fil) is the net polarization of the lattice 
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FIG. 1: (Color online) Schematic illustration of the unit cell 
of the cubic dipole lattice with embedded light source. 



dipoles and the emitter: 

P = d a S(r - r ) + ^2pjS(r - r 3 
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FIG. 2: (Color online) Isofrequency curves in xz plane, cal- 
culated for different dipole polarizabilities ao, zz - Normalized 
polarizability / a 3 is indicated near each curve. Calcu- 

lated was performed at qa — 0.15-7T. 



All the dipoles Pj are characterized by the identical po- 
larizability tensor a 



Pj = aE eXb (rj) , 



(3) 



Our goal is to determine the total electric field and 
polarizations, induced in the structure by the radiat- 
ing dipole do. This procedure includes field expansion 
over the Bloch eigenmodes with wavevectors fc, for which 
Eqs. ([l]),([2]) are independent. The results in a coordinate 
space are obtained by inverse Fourier transformation. In 
particular, the polarizations of lattice dipoles read 



Pj 



(BZ) 



(2tt) 3 



1-G(fc)a G , k (~r )d 



where Vq — a 3 is the unit cell volume, the integra- 
tion takes place over the Brillouin zone \k m \ < it /a, 
m — x,y, z, and 1 is 3 x 3 unity matrix. Radiating dipole 
position t*o enters Eq. Q and thus determines the effi- 
ciency of the lattice excitation. The quantity C in Eq. Q 
is the tensor interaction constant of the lattice, defined 
as 35 



C(fc) 



lim[Go,k(r)-G (r)] 
i — >o 



2iq 3 



(5) 



where Go.fc is the Green function of the photon with 
Bloch vector fc, 



G 0)fe (r) 
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and Go is the free photon Green function 



G (r)= [<z 2 +VV]i- 



(6) 



(7) 



The infinite lattice sums Q may be found either by 
Ewald summation [37 or by a Floquet-type summa- 
tion [35] ■ We have used the approach from Ref. [35, since 
it is preferential for fast evaluation of the integral Q. 
Electric field in the structure is given by the sum of the 
waves emitted by all the dipoles , 



E(r) = G (r - r )p 



)Pj 



(8) 



Eq. {[sj , by definition provides the Green function for the 
source embedded in the dipole lattice. Second term in 
Eq. ^ is given by Eq. Q where e lfcrj is replaced by 
Go, fe (r). 

From now we restrict ourselves to the case of uniaxial 
dipoles, when the only non-zero component of the tensor 
q is a zz . We assume that the the radiating dipole d 
is also directed along z axis. The TM-polarized Bloch 
eigenmodes of the system with given wavevector fc are 
found [S3 El] from the poles of Eq. Q 



1 



C(fc) = 0. 



(9) 



where G(fc) = C zz (k). Note, that Eq. ([9} is real for 
vanishing losses, because the imaginary part of the inter- 
action constant ^ cancels out with the radiative decay 
term in the polarizability: 



1 

OtQ,zz 



2iq 3 



1 



(10) 



Here «o,zz is the so-called bare dipole polarizability, cal- 
culated neglecting radiative decay [38] . Effective medium 
approximation for the solutions of Eq. ^ are the extraor- 
dinary TM-polarized modes, with the dispersion given 
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FIG. 3: (Color online) Spatial distribution of the dipole mo- 
ments \pz(r)\/po in the elliptic regime with Qo, 22 = a 3 /(47r). 
Insets schematically illustrate the geometry and the isofre- 
quency surfaces in wavevector space. Calculation was per- 
formed at qa = 0.15-7T and ro = 0.5az. 



by 



7 = + fc 



(11) 



Here e zz is the Maxwell-Garnett effective dielectric con- 
stant of the hyperbolic medium 



S zz = 1 



1 



(12) 



V/{A^ ZZ ) - 1/3 ' 
in the same approximation e xx — e yy = 1. 

III. DISPERSION AND GREEN FUNCTION 




FIG. 4: (Color online) Spatial distribution of the dipole mo- 
ments \pz(r)\/po in the hyperbolic regime, excited by the 
point emitter. Panels (a) and (b) show the distribution in 
the planes y = and x — y, respectively. Insets schemat- 
ically illustrate the geometry and the isofrequency surfaces 
in wavevector space. Calculation was performed at «o,zz = 
— 6a 3 /(47r) and the same other parameters as Fig. [3]. 



in the effective medium model when nonlocal effects are 
taken into account [4TJ1 HTj . 



In this Section we first discuss the dispersion of the 
Bloch waves in the dipole lattice (Sec. Ill A I and then 
investigate in detail the emission pattern of the dipole 
embedded in the lattice (Sec. IIIB). 



A. Isofrequency curves 

Isofrequency curves in the (k z ,k x ) plane, found from 
Eq. (J9J) for different polarizabilities ao,zz, are shown 
on Fig. [5] Depending on the polarizability, dispersion 
curves are either elliptic or hyperbolic, in agreement 
with Eq. ( fl2] ). The curves are generally well described 
by the effective medium approximation (11). However, 
an intermediate "mixed" regime is possible for ao.zz ~ 
— 1.3a 3 /(47r) (blue dashed curve), when two Bloch modes 
with hyperbolic and elliptic dispersion coexist in the 
structure. Such isofrequency curves can not be described 
by the Maxwell-Garnett theory Eqs. ( 11 ),( 12 1, which pre- 



dicts only one TM mode for given frequency. They were 
analyzed in Ref. |35] in more details and can be obtained 



B. Green function 

Here we will focus on the spatial distribution Q of 
the dipole moments |p(tj-)| in the discrete lattice un- 
der the point dipole excitation. Results of calcula- 
tion for the dipole polarizabilities «o,zz = a 3 /(47r) and 
&o,zz — — 6a 3 /(47r), corresponding to elliptic and hyper- 
bolic regimes, are shown on Fig. |] and Fig. gj respec- 
tively. Calculation demonstrates that the distribution is 
qualitatively different in hyperbolic regime: the pattern 
is strongly anisotropic and has characteristic cross-like 
shape, see Fig. |4j Moreover, in hyperbolic case the pat- 
tern depends on the azimuthal direction: it has distinct 
vertical ripples in the plane y — (Fig. [3^,), which are 
absent in the plane y = x (Fig. |3|d) . 

To understand these results it is instructive to com- 
pare them with Green function in the effective medium 
approximation, see Fig. [5j This approximation allows 
one to obtain the solution in a closed form [51 132]. In 
the case of a vertical orientation of the radiating dipole, 
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FIG. 5: (Color online) Spatial distribution of the polarization 
PeB,z(r)/(poa 3 ) induce by the point source in (a) effective 
elliptic medium with e zz = 2.5 and (b) effective hyperbolic 
medium with e zz = —1. Insets schematically illustrate the 
geometry and the isofrequency surfaces in wavevector space. 



xz plane. 
Fig. |3| 



Other calculation parameters are the same as for 



Po 



Green function reads 



E eS (r) = (g 2 + VV) ft £- 



jqFt. 



R = V £ zz(x - x ) 2 + e zz (y - y ) 2 + (z - z ) 2 . (13) 

This is generalization of Eq. ^ in the case of uniaxial 
medium. The relation between electric field and polar- 
ization in the effective medium model is local, 



47rP cff = (e cff - l)E cS 



It should be stressed that the effective medium approx- 
imation is not applicable on the spatial scales smaller 
than the lattice constant a. Consequently, the problem 
of point radiating dipole in discrete structure can not 
be reduced to the effective medium one. The effects of 
the radiating dipole position within the unit cell are also 
beyond the effective medium approximation. Thus, the 
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FIG. 6: Isofrequency curves of the dipole lattice in the hyper- 
bolic regime. Solid and thin lines correspond to numerical cal- 
culation and effective medium approximation Eq. (111. Other 
parameters are the same as for Fig. [4] The inset schemati- 
cally indicates the Brillouin zone of the square lattice, point 
E corresponds to k x = k y = 7r/(v2«) • 



results of two models, discrete and effective, may be com- 
pared only qualitatively. 

In order to clarify the ambiguity we have evaluated on 
Fig. [5] the polarization (14) at the discrete set of square 



lattice points rj and the radiating dipole is located at 
the point Tq — 0.5az, see the inset of panel (a). Fig. [5^, 
and Fig. [Hb show the spatial distribution of the polariza- 
tion Eq. ( 14 ) for the values effective dielectric constants 
off e zz = 2.5 and e zz = —1, corresponding to Fig. [3] and 
Fig. [4] In the case of the material with elliptic disper- 
sion the emission pattern is qualitatively the same as for 
the isotropic medium. The near field is concentrated at 
the emitter origin, r = ro, while the far-field is emitted 
perpendicularly to the dipole axis. The pattern changes 
dramatically in the hyperbolic case (Fig. [5b). The dis- 
tribution has a distinct cross-like shape, typical for hy- 
perbolic medium [5J [TT]. In the elliptic case, the only 
field singularity is that at the origin R = \r — Vq\ = 0. 
In the hyperbolic medium this singular point becomes a 
conical surface, where the field intensity is concentrated. 
Radiated waves are propagating within the cone R 2 > 0, 
and are evanescent outside this cone. Energy flow di- 
rections are normal to the isofrequency surfaces, so such 
cone in r-space is a direct counterpart of the hyperbolic 
dispersion curves in fc-space. 

Comparing numerical and effective medium results, 
Fig. [3] and Fig. [3^,, we see that in the elliptic case the 
Green function is qualitatively the same as in the effec- 
tive medium approximation. Weak spatial modulation of 
(14) the dipole moments |p(rj)|, seen on Fig. [3) is related to 



the deviations from the effective medium theory Eq. ( 13 1, 
which, as mentioned above, is not completely valid for 
the point excitation. Distinct cross-like distribution of 
Fig. [4^, is a fingerprint of the hyperbolic regime, similar 
to the effective medium approximation of Fig. [5}d. Com- 
paring Fig. [5b and Fig. [iK, we see, that in the discrete 



case the singularity in the effective medium solution ( 13 ) 
at the conical surface R — is smeared out and even 
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FIG. 7: (Color online) Spatial distribution of the dipole mo- 
ments \p z (r)\/po in the mixed hyperbolic-elliptic regime with 
C-o,zz = — 1.3a 3 /(47r). Insets schematically illustrate the ge- 
ometry and the isofrequency surfaces in wavevector space. 
Other calculation parameters are the same as for Fig. [5] 



vanishes at large enough distances, where the effective 
approximation also breaks down. This is qualitatively 
explained by the presence of the wavevector cutoff ~ n fa- 
Tine second striking difference between Fig. [4] and its ef- 
fective medium counterpart Fig. [5]d is the strong spatial 
modulation of the distribution in the y = plane, man- 
ifested as vertical ripples. Such modulation is obviously 
beyond the effective medium approximation of Fig. [5]d. 
In particular, the ripples on Fig. |4|a) turn out to be the 
result of the interference of the Bloch waves with wavec- 
tors k x = ±7r/o, corresponding to the boundaries of the 
Brillouin zone. To check this hypothesis we have plot- 
ted on Fig. [6] the isofrequency curves in T— X and T—M 
directions. Interference pattern in the planes y = and 
y = x should depend on the dispersion along T— X and 
T—M, respectively. Since dk z /dk x = at k x = ir/a (right 
panel of Fig. [6]) , there is a singularity in the density of 
states propagating along x direction, promoting the ver- 
tical ripples. This singularity is absent for the T—M di- 
rection, where the behavior of the isofrequency curves 
at the Brillouin zone edge is different. Panels (a) and 
(b) of Fig. [i] present the spatial distribution Q of the 
dipole moments [p(r,-)| in the planes y = and y = x, 
respectively. We see that the spatial modulation in the 
plane y = x is absent, cf. Fig. [4^, and Fig. |4Jd, which 
corroborates our explanation. The discovered effect can 
be thought of as the manifestation of the Van Hove band 
edge singularity [33] in the Green function [4"4] . 

Discrete Green function calculated for the dipole po- 
larizability ao. zz = — 1.3a 3 /(47r), corresponding to mixed 
elliptic- hyperbolic regime, is shown on Fig.[7| In this case 
the ripples are absent, because the isofrequency curves do 
not reach the Brillouin zone boundary k x — ±w/a, see 
dashed curve on Fig. [2] Far-field emission along x direc- 
tion is possible due to the modes with elliptic dispersion, 
providing weak background to the field of the hyperbolic 
modes. 



IV. PURCELL FACTOR 

Here we investigate the role of the discreetness on the 
Purcell factor determining the characteristics of the spon- 
taneous emission of the radiating dipole inside the ma- 
terial. The Purcell factor / and the Lamb shift / for 
the radiating dipole can be found from the Green func- 
tion ^ , evaluated at the dipole origin [H |SJ Hlj , see also 

Ref.aa 



/ + u = 1 



3i£z(r ) 
2q 3 Po 



(15) 



Here the dimensionless Lamb shift I is formally under- 
stood as a radiative correction to the resonance frequency 
of the radiating two-level system, normalized to its free 
space decay rate. Gathering Eqs. ( 15 ) , ([8]) , Q together, 
we find the result in a compact form 



f + il = 



2q 3 J (2tt)3 l/a zz ~ C(k) - 10 

(BZ) 



(16) 



The frequency w, entering the wavevector q in Eq. (16 1, 



is determined by the transition energy of the radiating 
dipole do. It is clear from the structure of Eq. (16), that 



the Purcell factor is determined by the pole contribution, 
corresponding to emission of photons with the dispersion 
given by Eq. ([9|. We note, that the first term in right- 
hand-side of Eq. (15) has canceled out in Eq. (16) with 



the pole contribution in the free space Green function 
Gk,zz(ro) at q = k. We stress, that despite the classi- 
cal formulation of the problem, the results for the emis- 
sion rate and photon Green function may be equivalently 
obtained by the quantum-mechanical calculation, either 
using the Fermi Golden rule [5] or the local quantization 
framework [4]. 

Numerical results for the dependence of the Purcell 
factor on the radiating dipole position within the unit 
cell of the structure are presented in Figs. |HJ Fig. [9] 
Fig. [TO] shows the frequency dependence of the Purcell 
factor. Figures demonstrate that the Purcell factor is 
much larger in hyperbolic regime than in the elliptic one. 
It is very sensitive to the dipole position and strongly 
increases when the dipole approaches the lattice nodes. 
Before discussing these results in more details it is in- 
structive to compared them with the analytical theory. 



A. Analytical results 

Here we focus on the Purcell factor in the quasi-static 



limit q <C 7r/a. Eq. (16) can be then reduced to 



2g 3 



(27T) 



x u,n>y 



dC{k) 



dk z 



(17) 
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FIG. 8: (Color online) Purcell factor in the (a) hyperbolic 
and (b) elliptic regime as function of the source coordinate 
zq for xo — yo = 0. Thick solid black, thin solid red, and 
dashed black curves correspond to numerical calculations, the 



analytical results of Eq. (20 1 (panel a), Eq. (23 1 (panel b), and 



to a single dipole with corresponding polarizability (Eq. ( 25 1) 



respectively. Dashed curves are results (17 1 for source near 



single dipole at r = 0. Other parameters as the same as for 

Kg. a 



where the interaction constant in the effective medium 
approximation reads |36j 



G(fc) 



47r q 2 



V k 2 - q 2 



47T 

w 



2iq 3 
~3~ 



(18) 



The integral over k z in ( 16 ) is determined by the residues 
at the wavevectors zLk z (k x , k y ), being the solutions of 
Eq. (9) at given frequency. Integration over k x and k y 
in Eq. (17 1 is restricted to those vectors within the two- 
dimensional Brillouin zone, for which such solution ex- 
ists. The quantity G zz . sta t( r o) m ( jlTP is the quasistatic 
approximation of the Green function (JsJ) : G ZZ)S t a t (To) = 
Gk,zz(ro)\k=o,q=o- The value of G ZZjSt at is determined by 



(a) Infinite lattice 



(b) Single dipole 




FIG. 9: (Color online) (a) Purcell factor in hyperbolic medium 
as function of the source position in the unit cell, (b) Calcula- 
tion in single-dipole approximation Eq. ( 25 1. Calculation was 



carried out at j/o = and the same other parameters as for 
Fig. [4] Radiating dipole coordinates change within the square 
< xo < 1, < Zo < 1. Colors correspond to logarithmic 
scale, identical for both panels. 




FIG. 10: (Color online) Purcell factor in (a) hyperbolic and 
(b) elliptic regime as function of the frequency qa for ro = 
0.25ai. Notation and other parameters as the same as for 
Fig. [8] 



the near field of the lattice dipoles, closest to the radiat- 
ing one. Maximum Purcell factor can be then expected 
when the source is located on the vertical edge of the 
elementary cell, i.e. xq = yo = 0, zq ^ 0. In this case 
G zz , s tat(zo) can be reduced to 



G 



zz, st&t 



Oo) ~ zs + 



z a ( a - z o) 3 



(19) 



and grows dramatically when the emitter approaches the 
lattice node. Evaluating the derivative in Eq. (17) by 



means of Eq. ( 18 1 and performing the integration, we 



obtain the analytical result for the Purcell factor 



/hyp — 



32vr 2 



\V a G zz , stat (z)f 



(20) 



in the hyperbolic medium. Here /c Zjmax 3> q is the cutoff 
for the wavevector k z , existing due to the finite extent 
of the Brillouin zone. Its value depends on the effective 
dielectric constant. 




1 < e < 



< -1 . 



(21) 



Thus, Eq. (20) provides a compact analytical result for 



the Purcell factor in the lossless hyperbolic medium. Its 
general structure can be understood as follows: the factor 
(&z,max/<?) 3 ~ l/(9a) 3 describes the enhancement of the 
photonic density of states as compared to the vacuum. 
The second factor |VoG ZZjSt at(-z)| 2 reflects the coordinate 
dependence of the Purcell factor, governed by the near- 
field of the neighboring dipoles. Near the lattice nodes 
Eq. ( 20 ) can be simplified to 



Ayp(g,Z-X))» 8g3| ^ |6 



(22) 
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where we assumed that |e| < 1. 

Similar calculation can be also performed in the ellip- 
tic case, when e zz > 0. It should be noted, that in the 
effective medium approximation, the Purcell factor for 
the axial dipolc in the elliptic medium is unity, indepen- 
dently of the value of e zz [21] • Local-field corrections can 
still promote high decay rate. The answer reads 



/ell = \G S 



V (s - 1) 



4tt 



(23) 



This expression depends on the local field intensity, sim- 



ilarly Eq. ( 20 1 , but is smaller by the factor 



/hyp 

/oil 



k 3 



2q 3 



(24) 



since the density of states in elliptic medium is smaller. 
In order to distinguish between the local field effects and 
the collective effects due to density of states enhancement 
in the medium it instructive to analyze also the Purcell 
factor for a source located in vacuum near single dipole 
in the point r = 0. The result reads [37] 



f 1 = l + ^lm[a zz Gl zz (r )], (25) 
here the second term is the field of the emitter, reflected 



from the dipole. In the quasistatic limit q — > Eq. ( 25 ) 
reduces to 



1 + 



|*3| 



(26) 



Both Eqs. (22 1 and (26) demonstrate divergency when 
z tends to zero. However, their dependence on the 
wavevector q is quite different: Eq. (22) diverges as 
1/q 3 at small q, while Eq. (26) does not depend on q 
at all. This divergency is a characteristic effect of pho- 
tonic density of states enhancement in the hyperbolic 
medium [HI [Ml US]- 



B. Numerical results 

Now we discuss the calculated dependence of the Pur- 
cell factor on the source position and on the transition 
frequency u> = cq, shown on Fig. [8] Fig. [9] and Fig. [TO] 
The calculation confirms the singular behavior of the 
Purcell factor in the hyperbolic case when the source ap- 



i). The 



proaches the lattice nodes (solid curve on Fig. 
singularity is excellently described by analytica 
(thin red curve). The interaction of the emitter with 



Eq. (20) 



a single dipole (Eq. (25)) provides substantially smaller 
enhancements (dashed black curve), although it also di- 
verges at z — and z = a. Additional comparison of 
the exact calculation in medium with the result for a sin- 
gle dipole is presented by the Purcell factor dependence 
on the two coordinates x and z in the unit cell, shown 




FIG. 11: (Color online) Purcell factor (solid lines) and Lamb 
shift (dashed lines) dependence on the coordinate zo of the 
source in the unit cell for different values of qa. Calculation 
parameters as the same as for Fig. [4] 



on Fig. [9} Eq. (25), taking into account only single lat- 
tice dipole at the point r — 0, satisfactory reproduces 
the Purcell factor pattern near this point. Corresponding 
two-dimensional plot of the Purcell factor in the quadrant 
< x, y < 1 is shown in Fig.[9^b). Comparing two panels 
of Fig. [9] we see, that near the corner r = the angu- 
lar dependence is approximately given by (3 cos 2 9 — l) 2 , 
where 9 is a polar angle. Still, single dipole model with 
corresponding polarizability ao. zz = — 6a 3 /(47r) consid- 
erably underestimates the absolute values of the Pur- 
cell factor. The satiation is different in the elliptic case, 
where all three approaches, namely numerical calcula- 
tion according to Eq. (16), single dipole model (25) with 
&o,zz = a 3 /(47r) and analytical model (23) provide simi- 
lar results, see Fig. [8]d. 

The failure of the single dipole model in the hyperbolic 
medium is also revealed in the frequency dependence of 



the Purcell factor, Fig. 10 1. The dashed curve, calculated 
for single dipole, tends to the limit (26), which is fre- 



quency independent. However, the Purcell factor in the 
hyperbolic medium diverges at low frequencies as 1/q 3 , 
according to Eq. (22). The Lamb shift I, calculated in 



hyperbolic medium for different values of qa is presented 
in Fig. [IT] by the dashed curves. Lamb shift is of the 
same order as the Purcell factor (dashed curve) and has 
similar near-field singularities at the node sites. 

To summarize, Figs. [8]- [ll| underline the importance 
of the local-field effects in the hyperbolic medium and 
confirm the collective origin of the spontaneous emission 
enhancement. 



V. CONCLUSIONS 

We have developed the analytical theory of light- 
matter coupling in discrete hyperbolic metamaterials in 
the framework of the discrete model of a cubic lattice 
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of uniaxial resonant dipoles. We have calculated Pur- 
cell factor, Lamb shift, and Green function for such a 
discrete model, and we have demonstrated that the opti- 
mal emitter position is in the local field maxima, close to 
the lattice nodes. We have demonstrated that the den- 
sity of states is drastically enhanced in the hyperbolic 
regime as compared to other cases including vacuum, el- 
liptic regime, or single resonant dipole case. As a result, 
a huge number of lattice dipoles are efficiently excited 
by the emitter, which has been visualized by calculating 
the Green function of the lattice. The Green function 
has a shape of a conus: the field propagates along the 
directions close to symmetry axis z and decays in the xy 
plane. Discrete character of the problem results in the 
strong spatial modulation of the Green function. 

The calculated absolute values of the Purcell factor are 
rather challenging for the current realization of metama- 
terials. This is mainly due to the point dipole approx- 
imation we have utilized: as distance to the scatterers 
becomes comparable to their sizes, higher order multi- 



poles must be accounted for. This will inevitably reduce 
the local field and suppress the Purcell factor. Achieving 
large density of states enhancement described by Eq. (p4|) 



is also not easy. Finally, the losses are inevitable and can 
significantly influence the numerical answers. Neverthe- 
less, we believe that our results will remain qualitatively 
correct for more complex settings, and they provide an 
important insight into the physics of hyperbolic metama- 
terials. 
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